
##########Marginal Means - Replication - Spilker et al. 2020

library(foreign)
library(cregg)
library(ggplot2)
library(grid)
library(gtable) 

###Analysis Vietnam
####General Marginal Means
vietnam <- read.dta("conjoint_data_R.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


###Look at province 1 (=Hanoi) only
vietnam <- read.dta("conjoint_data_R_procode1.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


###Look at province 74 - Thu Dau Mot only
vietnam <- read.dta("conjoint_data_R_procode74.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


#### Look at Province 79 - Ho Chi Minh City only
vietnam <- read.dta("conjoint_data_R_procode79.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)



#Low income

vietnam <- read.dta("conjoint_data_R_V_poor.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)



#High income
vietnam <- read.dta("conjoint_data_R_V_rich.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)



#education high
vietnam <- read.dta("conjoint_data_R_V_educ_high.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


#educ low
vietnam <- read.dta("conjoint_data_R_V_educ_low.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


#young
vietnam <- read.dta("conjoint_data_R_V_below_30.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


#between 30 and 50
vietnam <- read.dta("conjoint_data_R_V_between_30_50.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


#above 50
vietnam <- read.dta("conjoint_data_R_V_above_50.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


#Climate Change deniers

vietnam <- read.dta("conjoint_data_R_V_CC_deniers.dta") 

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


#Climate Change believers

vietnam <- read.dta("conjoint_data_R_V_CC_believers.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(vietnam, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)






###Kenya
kenya <- read.dta("conjoint_data_R_Kenya.dta")


m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)



###Kenya: ethnicity 1 - Kamba

kenya <- read.dta("conjoint_data_R_Kenya_Kamba.dta")


m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


###Kenya: ethnicity 2 - Luo

kenya <- read.dta("conjoint_data_R_Kenya_Luo.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


###Kenya: ethnicity 3 - Luhya

kenya <- read.dta("conjoint_data_R_Kenya_Luhya.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


###Kenya: ethnicity 4 - Kikuyu
kenya <- read.dta("conjoint_data_R_Kenya_Kikuyu.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)



###Kenya: Nairobi

kenya <- read.dta("conjoint_data_R_Kenya_Nairobi.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)



###Kenya: Mombasa
kenya <- read.dta("conjoint_data_R_Kenya_Mombasa.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)



###Kenya: Kisumu
kenya <- read.dta("conjoint_data_R_Kenya_Kisumu.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)



###Kenya: Low income
kenya <- read.dta("conjoint_data_R_Kenya_poor.dta")


m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


###Kenya: High income
kenya <- read.dta("conjoint_data_R_Kenya_rich.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)

###Kenya: Highly educated
kenya <- read.dta("conjoint_data_R_Kenya_educ_high.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)



###Kenya: Education low
kenya <- read.dta("conjoint_data_R_Kenya_educ_low.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)



###Kenya: Climate change believers
kenya <- read.dta("conjoint_data_R_Kenya_CC_believers.dta")


m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)

###Kenya: Climate change deniers
kenya <- read.dta("conjoint_data_R_Kenya_CC_deniers.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)




###Kenya: Residents between 30-50
kenya <- read.dta("conjoint_data_R_Kenya_between_30_50.dta")


m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)

###Kenya: Residents above 50
kenya <- read.dta("conjoint_data_R_Kenya_above_50.dta")

m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)


###Kenya: Residents below 30

kenya <- read.dta("conjoint_data_R_Kenya_below_30.dta")


m1<-accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON
m2 <- cj(kenya, m1, id = ~ respondent_num, estimate = "mm", h0 = 0.5)
# plot MMs
plot(m2, vline = 0.5)

















## subgroup analysis
kenya$ethnicity <- cut(kenya$Q3, 4)
edit(kenya)
x <- cj(na.omit(kenya), accept ~ GENDER + AGE + EDUCATION + INCOME + ETHNICITY + REASON,
id = ~ respondent_num, estimate = "mm", h0 = 0.5, by = ~ ethnicity)
plot(x, group = "ethnicity", vline = 0.5)









data("immigration")

## subgroup analysis
immigration$ethnosplit <- cut(immigration$ethnocentrism, 2)
x <- cj(na.omit(immigration), ChosenImmigrant ~ Gender + Education + LanguageSkills,
id = ~ CaseID, estimate = "mm", h0 = 0.5, by = ~ ethnosplit)
plot(x, group = "ethnosplit", vline = 0.5)

f1 <- ChosenImmigrant ~ Gender + Education + LanguageSkills + CountryOfOrigin + Job + JobExperience + JobPlans + ReasonForApplication + PriorEntry


d1 <- cj(immigration, f1, id = ~ CaseID, estimate = "mm", h0 = 0.5)
# plot MMs
plot(d1, vline = 0.5)

amce(residents_vietnam, accept_num ~ female_num + age_num + education_num + income_num + ethnicity_num + reason_num, id= ~ respondent_num)


amce(residents_vietnam, accept_num ~ female_num,  id= ~ respondent_num)

amce(residents_vietnam, accept ~ age,  id= ~ respondent_num)




